Regulatory control of somitogenesis
18 August, 2022
Having matched RNA and ATAC-seq profiles provides an opportunity to characterise the regulatory mechanisms that control somitogenesis.
## expression data
meta.rna <- read.table(paste0(dir, "RNA-seq/data/metadata_RNAseq.tsv"),
stringsAsFactors = FALSE, header = TRUE)
meta.rna <- meta.rna[meta.rna$QC == 1 & meta.rna$wrongStage == 0,]
meta.rna$stage <- paste0("stage", meta.rna$stage)
meta.rna$somite <- factor(meta.rna$somite, levels=c("SIII","SII","SI"))
data.rna <- read.table(paste0(dir, "RNA-seq/data/geneCounts.NORM_batchCorrected_14PCs.tsv"),
check.names = FALSE)
data.rna <- data.rna[,c('gene', meta.rna$sample)]
universe <- data.rna$gene
names(universe) <- row.names(data.rna)
## accessibility data
meta.atac <- read.table(paste0(dir, "ATAC-seq/data/metadata_ATACseq_GQ.tsv"),
stringsAsFactors = FALSE, header = TRUE)
meta.atac <- meta.atac[meta.atac$QCpass==1,]
meta.atac$stage <- paste0("stage", meta.atac$stage)
meta.atac$somite <- factor(meta.atac$somite, levels=c("SIII","SII","SI"))
data.atac <- readRDS(paste0(dir, "ATAC-seq/results/04_peakCounts_csawMerged.NORM.batchCorrected_18PCs.Rds"))Peak to gene linking
One way to do this is to identify chromatin loci in the vicinity of DE genes that show coordinated activity patterns; these loci can be postulated as putative regulatory elements driving the changes in gene expression. We can exploit the variation across development, which usually involves large changes in activity, to find gene-peak pairs using correlation analysis. To determine whether a given correlation score is significantly larger than expected by chance, we compare correlation values of the same gene against peaks that are in other chromosomes, which do not contribute to cis regulatory mechanisms.
For this analysis, we use the method implemented in FigR
(with minor modifications), which uses chromVAR to
determine a set of background peaks matched for accessibility and GC
content to model the null distribution of correlation scores for each
gene-peak pair. Naturally, the analysis is restricted to the 43 samples
where we generated good-quality libraries for both the RNA and ATAC
datasets. We look for correlated peaks within 100kb of a DE gene.
There are 15734 gene-peak pairs that are significantly correlated. These involve nearly 6K DE genes, linked to ~13.5K peaks.
c(n_genes = length(unique(links$Gene)),
n_peaks = length(unique(links$PeakRanges)))## n_genes n_peaks
## 5909 13657Most gene-peak pairs show moderate correlation values (median=0.382547). However, a few pairs with very low correlation coefficients are deemed significant by the method. These correspond to cases where gene expression and accessibility are not associated (correlation values close to 0) but the set of background peaks used to compute significance are biased towards negative correlations. This leads to a significant p-value even if the correlation of the pair in question is not significant itself.
ggplot(links, aes(1, rObs)) +
geom_violin() +
geom_boxplot(width=0.1) +
xlab("") +
ylab("Spearman correlation") +
th + theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
Thus, to consider a gene-peak pair as significantly correlated we require a minimum correlation coefficient of 0.3 on top of a significant p-value, taking the number of links to 12.8K.
links <- links[links$rObs>0.3,]
nrow(links)## [1] 12803Although most genes are linked to a single or a few peaks, a small proportion are connected to more than 5 peaks (349 genes). In contrast, the vast majority of peaks (88.91%) regulate a single gene.
par(mfrow=c(1,2))
plot(table(table(links$Gene)), bty="l",
xlab="number of peaks per gene",
ylab="number of genes")
plot(table(table(links$PeakRanges)), bty="l",
xlab="number of genes per peak",
ylab="number of peaks")
# summary(as.numeric(table(links$Gene)))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 1.00 1.00 2.00 2.452 3.000 30.000
# summary(as.numeric(table(links$PeakRanges)))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 1.000 1.000 1.000 1.147 1.000 8.000 The set of defined pairs involve around half of all the DE genes, for both those that are differential across somite trios, and across stages.
## check how many of the correlated regions are differential
dars.somite <- read.table(paste0(dir, "RNA+ATAC/results/01_DAregions_summary_somiteTrios_withClusters.tsv"))
dars.stage <- read.table(paste0(dir, "RNA+ATAC/results/02_DAregions_summary_stage_fate.tsv"))
dars <- union(row.names(dars.somite), row.names(dars.stage))
links$DE <- ifelse(links$Gene %in% degs.somite$gene & links$Gene %in% degs.stage$gene, "both",
ifelse(links$Gene %in% degs.somite$gene, "somite",
ifelse(links$Gene %in% degs.stage$gene, "stage", "nonDE")))
links$DE_stage <- degs.stage[match(links$Gene, degs.stage$gene),]$fate
# the same proportion of DE genes in each fate are linked
links$DA <- ifelse(links$PeakRanges %in% row.names(dars.somite) & links$PeakRanges %in% row.names(dars.stage), "both",
ifelse(links$PeakRanges %in% row.names(dars.somite), "somite",
ifelse(links$PeakRanges %in% row.names(dars.stage), "stage", "nonDA")))
links$DA_stage <- dars.stage[match(links$PeakRanges, row.names(dars.stage)),]$fate
df <- data.frame(n_DE_linked = c(length(unique(links[links$DE %in% c("both", "somite"),]$Gene)),
length(unique(links[links$DE %in% c("both", "stage"),]$Gene))),
DE_pct = c(round(length(unique(links[links$DE %in% c("both", "somite"),]$Gene))/nrow(degs.somite)*100,2),
round(length(unique(links[links$DE %in% c("both", "stage"),]$Gene))/nrow(degs.stage)*100,2)),
n_DA_linked = c(length(unique(links[links$DA %in% c("both", "somite"),]$PeakRanges)),
length(unique(links[links$DA %in% c("both", "stage"),]$PeakRanges))),
DA_pct = c(round(length(unique(links[links$DA %in% c("both", "somite"),]$PeakRanges))/nrow(dars.somite)*100,2),
round(length(unique(links[links$DA %in% c("both", "stage"),]$PeakRanges))/nrow(dars.stage)*100,2)))
row.names(df) <- c("somite", "stage")
dfAlthough some peaks are nearby their putative target, the majority are dozens of kilobases away.
## TSS coordinates
TSSg <- FigR::mm10TSSRanges
# distance between gene and peak
tmp <- unlist(lapply(strsplit(links$PeakRanges, ":"), '[[', 2))
links.gr <- GenomicInteractions(GRanges(seqnames = unlist(lapply(strsplit(links$PeakRanges, ":"), '[[', 1)),
IRanges(start = as.numeric(unlist(lapply(strsplit(tmp, "-"), '[[', 1))),
end = as.numeric(unlist(lapply(strsplit(tmp, "-"), '[[', 2))))),
TSSg[match(links$Gene, TSSg$gene_name),])
distance <- distance(anchorOne(links.gr), anchorTwo(links.gr))
links$distance <- distance
bins <- data.frame(distance=c('0-200','200-2e3','2e3-5e3','5e3-10e3','10e3-50e3','50e3-100e3'),
n=c(sum(distance <= 200),
sum(distance > 200 & distance <= 2e3),
sum(distance > 2e3 & distance <= 5e3),
sum(distance > 5e3 & distance <= 10e3),
sum(distance > 10e3 & distance <= 50e3),
sum(distance > 50e3 & distance <= 100e3)))
bins$n_pct <- round(bins$n/length(distance)*100,2)
binsThe median distance between linked peaks and genes is of 37.845kb.
Supporting the regulatory potential of the linked peaks, a larger proportion overlap ENCODE cCREs characterised as enhancers (H3K27ac-high, H3K4me3-low) compared to all peaks detected in somites, and all peaks within 100kb of a DE gene. Furthermore, the overlap greatly increases if we restrict to links with the strongest correlations.
encode <- read.table(paste0(dir, "ATAC-seq/results/05_peaks_overlapExternalData.tsv"), sep="\t")
row.names(encode) <- paste0("chr", row.names(encode))
encode$encode <- ifelse(encode$pELS | encode$dELS, "enhancer",
ifelse(encode$PLS, "promoter",
ifelse(encode$CTCF_bound, "CTCFonly",
ifelse(encode$H3K4me3, "H3K4me3only", "not_in_encode"))))
links$encode <- encode[match(links$PeakRanges, row.names(encode)),]$encode
links$fantom <- encode[match(links$PeakRanges, row.names(encode)),]$FANTOM
links$vista <- encode[match(links$PeakRanges, row.names(encode)),]$VISTA
links$encode <- as.factor(links$encode)
# round(rbind(all_peaks = prop.table(table(encode$encode)[1:4])*100,
# peaks_proximal_to_DEGs = prop.table(table(encode[unique(cisCor$PeakRanges),]$encode)[1:4])*100,
# peaks_linked_to_DEGs = prop.table(table(links$encode)[1:4])*100,
# linked_0.5 = prop.table(table(links[links$rObs>0.5,]$encode)[1:4])*100,
# linked_0.65 = prop.table(table(links[links$rObs>0.65,]$encode)[1:4])*100), 2)[,c(4,2,1,3)]
df <- data.frame(type = rep(c("enhancer", "not_in_encode"), each=7),
prop = c(sum(encode$encode=="enhancer")/nrow(encode)*100,
sum(encode[unique(cisCor$PeakRanges),]$encode=="enhancer")/nrow(encode[unique(cisCor$PeakRanges),])*100,
sum(links[links$rObs>0.3,]$encode=="enhancer")/nrow(links[links$rObs>0.3,])*100,
sum(links[links$rObs>0.4,]$encode=="enhancer")/nrow(links[links$rObs>0.4,])*100,
sum(links[links$rObs>0.5,]$encode=="enhancer")/nrow(links[links$rObs>0.5,])*100,
sum(links[links$rObs>0.6,]$encode=="enhancer")/nrow(links[links$rObs>0.6,])*100,
sum(links[links$rObs>0.7,]$encode=="enhancer")/nrow(links[links$rObs>0.7,])*100,
sum(encode$encode=="not_in_encode")/nrow(encode)*100,
sum(encode[unique(cisCor$PeakRanges),]$encode=="not_in_encode")/nrow(encode[unique(cisCor$PeakRanges),])*100,
sum(links[links$rObs>0.3,]$encode=="not_in_encode")/nrow(links[links$rObs>0.3,])*100,
sum(links[links$rObs>0.4,]$encode=="not_in_encode")/nrow(links[links$rObs>0.4,])*100,
sum(links[links$rObs>0.5,]$encode=="not_in_encode")/nrow(links[links$rObs>0.5,])*100,
sum(links[links$rObs>0.6,]$encode=="not_in_encode")/nrow(links[links$rObs>0.6,])*100,
sum(links[links$rObs>0.7,]$encode=="not_in_encode")/nrow(links[links$rObs>0.7,])*100))
ggplot(df, aes(rep(1:7,2), prop, colour=type)) +
geom_point(size=2) +
geom_line() +
xlab("") +
ylab("% of peaks") +
ylim(c(0,100)) +
scale_x_continuous(breaks=1:7,
labels=c("all peaks", "peaks near DEGs", paste0("linked peaks > 0.", 3:7))) +
th + theme(axis.text.x = element_text(angle=45, hjust=1),
panel.grid.major.y = element_line(size=0.5, colour="grey90"),
panel.grid.minor.y = element_line(size=0.25, colour="grey95"))
Similarly, while only ~10% of all peaks overlap with FANTOM5 enhancers, the number increases to 14% for linked peaks and this goes up as the correlation threshold becomes more stringent.
df <- data.frame(type = "in_FANTOM",
prop = c(sum(encode$FANTOM)/nrow(encode)*100,
sum(encode[unique(cisCor$PeakRanges),]$FANTOM)/nrow(encode[unique(cisCor$PeakRanges),])*100,
sum(links[links$rObs>0.3,]$fantom)/nrow(links[links$rObs>0.3,])*100,
sum(links[links$rObs>0.4,]$fantom)/nrow(links[links$rObs>0.4,])*100,
sum(links[links$rObs>0.5,]$fantom)/nrow(links[links$rObs>0.5,])*100,
sum(links[links$rObs>0.6,]$fantom)/nrow(links[links$rObs>0.6,])*100,
sum(links[links$rObs>0.7,]$fantom)/nrow(links[links$rObs>0.7,])*100))
ggplot(df, aes(1:7, prop, colour=type)) +
geom_point(size=2) +
geom_line() +
xlab("") +
ylab("% of peaks") +
ylim(c(0,100)) +
scale_x_continuous(breaks=1:7,
labels=c("all peaks", "peaks near DEGs", paste0("linked peaks > 0.", 3:7))) +
th + theme(axis.text.x = element_text(angle=45, hjust=1),
panel.grid.major.y = element_line(size=0.5, colour="grey90"),
panel.grid.minor.y = element_line(size=0.25, colour="grey95"))
Highly-regulated genes
As mentioned above, although most genes are linked only to a single
or few peaks, a few hundred have many more connections. The authors of
FigR coined the term domains of
regulatory chromatin (DORCs) to refer to these loci of high
interactivity, and showed that they “are enriched for
lineage-determining genes and overlap with known super-enhancers”. In
our case, many of the genes linked to dozens of peaks correspond to late
Hox genes, which is expected since it has been shown that Hox expression
regulation relies on coordinated chromatin remodeling. Many other genes
are identified as well.
# Determine DORC genes
dorcGenes <- dorcJPlot(dorcTab = links,
cutoff = 6,
returnGeneList = TRUE)
This restricted set of genes with a large number of linked peaks are still enriched for key processes relevant for somitogenesis, such as patterning and the development and morphogenesis of the different somite derivative lineages. These data suggest that the expression of genes that are crucial in the development and differentiation of somites are under intricate regulatory control.
go <- topGOtable(DEgenes = dorcGenes,
BGgenes = universe,
topGO_method2 = "elim",
ontology = "BP",
geneID = "symbol",
fullNamesInRows = FALSE,
addGeneToTerms = TRUE,
mapping = "org.Mm.eg.db",
topTablerows = 100)
go[,c(2,4,7,9)]Interestingly, while genes most active at different fates are all as likely to be associated to at least one peak, there is a strong depletion of genes expressed in thoracic somites among DORC genes.
## check proportion per stage
tmp <- unique(links[,c('Gene','DE_stage')])
df <- rbind(all_DEGs = prop.table(table(degs.stage$fate))*100,
all_linked_DEGs = prop.table(table(tmp$DE_stage))*100,
DORC_DEGs = prop.table(table(tmp[tmp$Gene %in% dorcGenes,]$DE_stage))*100)
df[,c(1,4,2,3)]## cervical thoracic lumbar sacral
## all_DEGs 34.49304 17.902584 17.58449 30.01988
## all_linked_DEGs 32.70559 15.861059 17.99540 33.43796
## DORC_DEGs 38.05310 8.259587 15.63422 38.05310Regulatory control of DORC genes
Next, we can ask whether the peaks linked to the highly-regulated
genes are enriched for particular TF binding sites, compared to a set of
(matched) background peaks. These TFs have the potential of regulating
the expression of the gene in question. This hypothesis becomes stronger
if the expression of the TF itself is correlated with the accessibility
values of the peaks. FigR performs these two calculations,
and combines the results into a regulation score.
Positive/negative scores correspond to activating/repressing
interactions.
We observe a set of DORCs with very strong TFBS enrichment scores, all corresponding to enrichment for Nr6a1 binding sites; DORCs regulated by Nr6a1 show negative scores, indicating TF expression is negatively correlated with chromatin accessibility. Other transcription factors result in more modest motif enrichment scores, with both repressing and enhancing activities.
ggplot(fig.d, aes(Corr.log10P, Enrichment.log10P, color=Score, shape=as.factor(Motif=="Nr6a1"))) +
geom_point_rast(size=0.75) +
scale_color_gradientn(colours = jdb_palette("solar_extra"),
limits=c(-2,2),
oob = scales::squish,
breaks=scales::breaks_pretty(n=3)) +
geom_vline(xintercept = 0, lty=2, colour="grey") +
xlab(expression('(signed) -log'[10]*' p-value TF expression - accessibility correlation')) +
ylab(expression('-log'[10]*' p-value TF motif enrichment')) +
labs(shape = "Nr6a1 enrichment") +
th
Across all DORC genes, some TFs consistently show large regulation scores, as shown below, with Nr6a1 again standing out. Several of the TFs highlighted have been already picked up in previous analyses, given the enrichment of their binding sites in DA peaks.
rankDrivers(fig.d, rankBy = "meanScore", interactive = FALSE)
Generally, TFs are associated with a handful to a dozen different targets, but some factors show more widespread activity (again, Nr6a1, Ebf1 and Sox6, all acting as repressors, and several Hox genes).
rankDrivers(fig.d, score.cut = 1, rankBy = "nTargets", interactive = TRUE)## Ranking TFs by total number of associated DORCs ..## Using absolute score cut-off of: 1 ..## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## Please use `gather()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.We can then group DORCs (rows) that are regulated in a concerted manner by groups of TFs (columns).
score.cut <- 1.25
set.seed(2948)
regulons <- plotfigRHeatmap(figR.d = fig.d,
score.cut = score.cut,
row_names_gp = gpar(fontsize=8, fontface="italic"),
column_names_gp = gpar(fontsize=10),
show_row_dend = TRUE, k=4)
regulons <- draw(regulons)
## get cluster membership
DORCsToKeep <- fig.d %>% dplyr::filter(abs(Score) >= score.cut) %>% pull(DORC) %>% unique()
TFsToKeep <- fig.d %>% dplyr::filter(abs(Score) >= score.cut) %>% pull(Motif) %>% unique()
net.d <- fig.d %>% dplyr::filter(DORC %in% DORCsToKeep & Motif %in% TFsToKeep) %>%
reshape2::dcast(DORC ~ Motif) %>%
tibble::column_to_rownames("DORC") %>% as.matrix()
clusters <- as.data.frame(unlist(row_order(regulons)))
clusters$cluster <- paste0("cluster", substr(row.names(clusters), 1, 1))
clusters$gene <- row.names(net.d[clusters[,1],])
clusters$cluster <- factor(clusters$cluster, levels=paste0("cluster", names(row_order(regulons))))And visualise the expression of the DORC genes, to link each module to its activity pattern across stages.
data <- data.rna[match(clusters$gene, row.names(data.rna)),]
data <- t(apply(data, 1, function(x) (x-mean(x))/sd(x) ))
m <- meta.rna[match(colnames(data), meta.rna$sample),]
m$stage <- factor(m$stage, levels=paste0("stage", c(8,18,21,25,27,35)))
m$somite <- factor(m$somite, levels=c("SIII","SII","SI"))
order <- order(m$stage, m$somite)
## heatmap annotation
ha <- HeatmapAnnotation(df = data.frame(stage = m[order,]$stage,
somite = m[order,]$somite),
col = list(stage = cols.stage, somite = cols.somite))
Heatmap(data[,order],
cluster_columns = FALSE,
cluster_rows = FALSE,
col=colorRamp2(breaks = c(-4,-2,0,2,3),
colors = c("steelblue4","steelblue","white",
"indianred3","indianred4")),
row_names_gp = gpar(fontsize=8, fontface="italic"),
column_names_gp = gpar(fontsize=10),
name = "z-score",
show_row_names = TRUE,
show_column_names = FALSE,
top_annotation = ha,
row_split = clusters$cluster,
cluster_row_slices = FALSE)
Two main modules are evident, with genes expressed either early or late
in development. The role of Nr6a1 is prominent, and Hox genes
also show regulatory effects consistent with developmental
progression.
saveRDS(cisCor, paste0(dir, "RNA+ATAC/results/03_gene_peak_correlations.Rds"))
write.table(links, paste0(dir, "RNA+ATAC/results/03_gene_peak_links.tsv"), quote = FALSE, sep = "\t", row.names = FALSE)
saveRDS(fig.d, paste0(dir, "RNA+ATAC/results/03_FigR_network.Rds"))
saveRDS(net.d, paste0(dir, "RNA+ATAC/results/03_regulationScores_thr1.25.Rds"))
write.table(clusters, paste0(dir, "RNA+ATAC/results/03_regulons_clusters.tsv"), quote = FALSE, sep="\t", row.names = FALSE)sessionInfo()## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] BSgenome.Mmusculus.UCSC.mm10_1.4.0 BSgenome_1.60.0
## [3] rtracklayer_1.52.1 Biostrings_2.60.2
## [5] XVector_0.32.0 org.Mm.eg.db_3.13.0
## [7] FigR_0.1.0 motifmatchr_1.14.0
## [9] chromVAR_1.14.0 cowplot_1.1.1
## [11] Matrix_1.4-1 pbmcapply_1.5.0
## [13] doParallel_1.0.17 iterators_1.0.14
## [15] foreach_1.5.2 dynamicTreeCut_1.63-1
## [17] pcaExplorer_2.18.0 topGO_2.44.0
## [19] SparseM_1.81 GO.db_3.13.0
## [21] AnnotationDbi_1.54.1 graph_1.70.0
## [23] networkD3_0.4 circlize_0.4.14
## [25] ComplexHeatmap_2.8.0 RColorBrewer_1.1-3
## [27] BuenColors_0.5.6 MASS_7.3-56
## [29] GGally_2.1.2 ggrastr_1.0.1
## [31] ggpubr_0.4.0 ggrepel_0.9.1
## [33] ggplot2_3.3.5 dplyr_1.0.8
## [35] GenomicInteractions_1.26.0 InteractionSet_1.20.0
## [37] SummarizedExperiment_1.22.0 Biobase_2.52.0
## [39] GenomicRanges_1.44.0 GenomeInfoDb_1.28.4
## [41] IRanges_2.26.0 S4Vectors_0.30.2
## [43] BiocGenerics_0.38.0 MatrixGenerics_1.4.3
## [45] matrixStats_0.62.0
##
## loaded via a namespace (and not attached):
## [1] rappdirs_0.3.3 AnnotationForge_1.34.1
## [3] R.methodsS3_1.8.1 nabor_0.5.0
## [5] pkgmaker_0.32.2 tidyr_1.2.0
## [7] bit64_4.0.5 knitr_1.38
## [9] R.utils_2.11.0 DelayedArray_0.18.0
## [11] data.table_1.14.2 rpart_4.1.16
## [13] KEGGREST_1.32.0 TFBSTools_1.30.0
## [15] RCurl_1.98-1.6 AnnotationFilter_1.16.0
## [17] generics_0.1.2 snow_0.4-4
## [19] GenomicFeatures_1.44.2 RSQLite_2.2.12
## [21] tzdb_0.3.0 bit_4.0.4
## [23] webshot_0.5.3 xml2_1.3.3
## [25] httpuv_1.6.5 assertthat_0.2.1
## [27] DirichletMultinomial_1.34.0 viridis_0.6.2
## [29] xfun_0.30 hms_1.1.1
## [31] jquerylib_0.1.4 evaluate_0.15
## [33] promises_1.2.0.1 TSP_1.2-0
## [35] fansi_1.0.3 restfulr_0.0.13
## [37] progress_1.2.2 caTools_1.18.2
## [39] dendextend_1.15.2 dbplyr_2.1.1
## [41] Rgraphviz_2.36.0 igraph_1.3.0
## [43] DBI_1.1.2 geneplotter_1.70.0
## [45] htmlwidgets_1.5.4 reshape_0.8.9
## [47] Rmpfr_0.8-7 purrr_0.3.4
## [49] ellipsis_0.3.2 crosstalk_1.2.0
## [51] backports_1.4.1 annotate_1.70.0
## [53] gridBase_0.4-7 biomaRt_2.48.3
## [55] vctrs_0.4.1 ensembldb_2.16.4
## [57] Cairo_1.5-15 abind_1.4-5
## [59] cachem_1.0.6 withr_2.5.0
## [61] Gviz_1.36.2 checkmate_2.0.0
## [63] GenomicAlignments_1.28.0 prettyunits_1.1.1
## [65] cluster_2.1.3 seqLogo_1.58.0
## [67] lazyeval_0.2.2 crayon_1.5.1
## [69] genefilter_1.74.1 labeling_0.4.2
## [71] pkgconfig_2.0.3 vipor_0.4.5
## [73] ProtGenerics_1.24.0 seriation_1.3.5
## [75] nnet_7.3-17 rlang_1.0.2
## [77] lifecycle_1.0.1 miniUI_0.1.1.1
## [79] registry_0.5-1 filelock_1.0.2
## [81] BiocFileCache_2.0.0 doSNOW_1.0.20
## [83] GOstats_2.58.0 dichromat_2.0-0
## [85] rngtools_1.5.2 carData_3.0-5
## [87] base64enc_0.1-3 beeswarm_0.4.0
## [89] GlobalOptions_0.1.2 pheatmap_1.0.12
## [91] png_0.1-7 viridisLite_0.4.0
## [93] rjson_0.2.21 bitops_1.0-7
## [95] shinydashboard_0.7.2 R.oo_1.24.0
## [97] blob_1.2.3 shape_1.4.6
## [99] stringr_1.4.0 readr_2.1.2
## [101] jpeg_0.1-9 rstatix_0.7.0
## [103] shinyAce_0.4.1 ggsignif_0.6.3
## [105] CNEr_1.28.0 scales_1.2.0
## [107] memoise_2.0.1 GSEABase_1.54.0
## [109] magrittr_2.0.3 plyr_1.8.7
## [111] zlibbioc_1.38.0 threejs_0.3.3
## [113] compiler_4.1.0 BiocIO_1.2.0
## [115] clue_0.3-60 DESeq2_1.32.0
## [117] Rsamtools_2.8.0 cli_3.2.0
## [119] Category_2.58.0 htmlTable_2.4.0
## [121] Formula_1.2-4 tidyselect_1.1.2
## [123] stringi_1.7.6 shinyBS_0.61.1
## [125] highr_0.9 yaml_2.3.5
## [127] locfit_1.5-9.5 latticeExtra_0.6-29
## [129] sass_0.4.1 VariantAnnotation_1.38.0
## [131] tools_4.1.0 rstudioapi_0.13
## [133] TFMPvalue_0.0.8 foreign_0.8-82
## [135] gridExtra_2.3 farver_2.1.0
## [137] digest_0.6.29 FNN_1.1.3
## [139] pracma_2.3.8 shiny_1.7.1
## [141] Rcpp_1.0.8.3 car_3.0-12
## [143] broom_0.8.0 later_1.3.0
## [145] httr_1.4.2 biovizBase_1.40.0
## [147] colorspace_2.0-3 XML_3.99-0.9
## [149] splines_4.1.0 RBGL_1.68.0
## [151] plotly_4.10.0 gmp_0.6-5
## [153] xtable_1.8-4 poweRlaw_0.70.6
## [155] jsonlite_1.8.0 heatmaply_1.3.0
## [157] R6_2.5.1 Hmisc_4.7-0
## [159] pillar_1.7.0 htmltools_0.5.2
## [161] mime_0.12 NMF_0.24.0
## [163] glue_1.6.2 fastmap_1.1.0
## [165] DT_0.22 BiocParallel_1.26.2
## [167] codetools_0.2-18 utf8_1.2.2
## [169] lattice_0.20-45 bslib_0.3.1
## [171] tibble_3.1.6 curl_4.3.2
## [173] ggbeeswarm_0.6.0 gtools_3.9.2
## [175] magick_2.7.3 survival_3.3-1
## [177] limma_3.48.3 rmarkdown_2.13
## [179] munsell_0.5.0 GetoptLong_1.0.5
## [181] GenomeInfoDbData_1.2.6 reshape2_1.4.4
## [183] gtable_0.3.0